library(readxl)
library(tidyverse)
library(kableExtra)
library(sf)
library(mapview)
library(vegan)
library(devtools)
library(indicspecies)
library(SpadeR)
library(adespatial) #TBI, para comparar matrices de comunidad entre estaciones
library(RColorBrewer)
library(spdep)
download.file('https://raw.githubusercontent.com/biogeografia-master/scripts-de-analisis-BCI/master/biodata/funciones.R', 'funciones1.R', )
source('funciones1.R')
source('funciones2.R')
source('estadistica-zonal.R')
tabla_odk <- read_csv('tabla ODK.csv')
tabla_odk <- tabla_odk[grep('prueba', tabla_odk$codigo, ignore.case = T, invert = T), ]
tabla_odk$codigo[grep('h14m75', tabla_odk$codigo, ignore.case = T)] <- 'H10M75'
tabla_odk$codigo[grep('h1m14', tabla_odk$codigo, ignore.case = T)] <- 'H12M14'
tabla_odk$codigo <- toupper(tabla_odk$codigo)
tabla_odk$codigo[grep('h1m88', tabla_odk$codigo, ignore.case = T)] <- 'H2M88'
tabla_odk$codigo[grep('2h23h42', tabla_odk$codigo, ignore.case = T)] <- '2H23M42'
tabla_odk$codigo[grep('2h23h43', tabla_odk$codigo, ignore.case = T)] <- '2H23M43'
tabla_odk$codigo[grep('2H25M147', tabla_odk$codigo, ignore.case = T)] <- '2H24M147'
datos_orig <- read_xlsx(
path = 'Datos de muestreo Marzo y Mayo 2022.xlsx',
col_types = c('numeric', rep('text', 7), 'date', 'numeric', rep('text', 5)))
which(is.na(match(tabla_odk$codigo, datos_orig$Código))) #Todos los elementos tienen contrapartes entre tabla si "integer(0)"
## integer(0)
ncol(datos_orig) # Número de columnas
## [1] 15
nrow(datos_orig)
## [1] 255
colnames(datos_orig)
## [1] "Parcela" "Código" "Coordenadas" "Datum"
## [5] "Altura" "Incertidumbre" "Especie" "Colector/os"
## [9] "Fecha" "No. de individos" "Observaciones" "País"
## [13] "Provincia" "Municipio" "Área protegida"
datos_01 <- datos_orig %>%
mutate(lat = str_split(datos_orig$Coordenadas, pattern = ' ', simplify = T)[,1],
lon = str_split(datos_orig$Coordenadas, pattern = ' ', simplify = T)[,2]) %>%
mutate_at(c('lat', 'lon'), as.numeric)
datos_02 <- datos_01 %>%
select(-Parcela) %>%
mutate(Hexágono = as.numeric(gsub('.*H', '', gsub('M.*', '', Código)))) %>%
select(Hexágono, everything())
datos_02 %>%
select(Hexágono, Código, Especie) %>%
head() %>%
kable()
| Hexágono | Código | Especie |
|---|---|---|
| 13 | H13M1 | Pheidole subarmata |
| 13 | H13M2 | Dorymyrmex antillanus |
| 13 | H13M3 | Paratrechina longicornis |
| 13 | H13M4 | Pheidole subarmata |
| 13 | H13M5 | Solenopsis geminata |
| 13 | H13M6 | Pheidole jelskii |
datos_03 <- datos_02 %>%
inner_join(
tabla_odk %>%
select(codigo, `coordenadaid-Latitude`, `coordenadaid-Longitude`), by = c('Código' = 'codigo')) %>%
select(-lat,-lon) %>%
rename(lat = `coordenadaid-Latitude`,
lon = `coordenadaid-Longitude`)
datos_03[grep('^H13M1$', datos_03$Código, ignore.case = T), 'Hexágono'] <- 12 #Misplaced, corrected
datos_03$Código[grep('^H13M1$', datos_03$Código, ignore.case = T)] <- 'H12M1' #Misplaced, corrected
datos_03[grep('^H1M87$', datos_03$Código, ignore.case = T), 'Hexágono'] <- 4 #Misplaced, corrected
datos_03$Código[grep('^H1M87$', datos_03$Código, ignore.case = T)] <- 'H4M87' #Misplaced, corrected
datos_03$Especie <- gsub('Pheidole jelskii', 'Pheidole jelskii', datos_03$Especie)
datos_03$Especie <- gsub('Ontomachus ruginosa', 'Odontomachus ruginodis', datos_03$Especie)
datos_03$Especie <- gsub('Monomorium floricula', 'Monomorium floricola', datos_03$Especie)
datos_03$Especie <- gsub('Monomorium lanuginosa', 'Tetramorium lanuginosum', datos_03$Especie)
# datos_03 %>% write_csv('datos_03.csv')
# write_excel_csv(datos_03, 'datos_03_compatible_con_excel.csv')
# datos_03 <- read_csv('datos_03.csv')
datos_sf <- datos_03 %>%
st_as_sf(coords = c('lon', 'lat'))
st_crs(datos_sf) <- 4326
datos_sf %>% plot()
# st_write(datos_sf, 'datos_sf.gpkg', delete_dsn = T)
# Generar mapa de cuadros sin simbología
## Centroid
# datos_sf %>% mutate(i=1) %>%
# group_by(i) %>%
# summarize(geometry = st_union(geometry)) %>%
# st_centroid %>% mutate(lng )
## Mapa
mapa_cuadros <- mapView(
datos_sf,
col.regions = 'grey80',
alpha.regions = 0.3,
map.types = 'Esri.WorldImagery',
legend = F, zoom = 14,
zcol = 'Código') %>% leafem::addStaticLabels(label = datos_sf$Código) # %>%
# leaflet::setView(
# lng = -79.85136,
# lat = 9.15097,
# zoom = 15)
mapa_cuadros
# mapa_cuadros %>% mapshot(file = 'mapa_punto.png') #Genera archivo
# Capa de hexagonos
hex <- st_read('cuadricula-final.gpkg')
## Reading layer `cuadricula-final' from data source
## `C:\Users\User\Desktop\Universidad\Tesis\Datos de los muestreos\cuadricula-final.gpkg'
## using driver `GPKG'
## Simple feature collection with 28 features and 9 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 402697.2 ymin: 2041872 xmax: 403143.5 ymax: 2042260
## Projected CRS: WGS 84 / UTM zone 19N
hex %>% inner_join(#NUMERO DE NIDOS POR HEXAGONO
datos_03 %>%
group_by(Hexágono) %>%
summarise(n = n()), by = c('ENLACE' = 'Hexágono')) %>%
mapview(zcol = 'n')
datos_03 %>% group_by(Hexágono) %>%
summarise(n = n()) %>% arrange(desc(n)) %>%
head() %>%
kable()
| Hexágono | n |
|---|---|
| 12 | 19 |
| 16 | 16 |
| 21 | 16 |
| 20 | 15 |
| 13 | 13 |
| 6 | 12 |
# datos_03 <- read_csv('datos_03.csv') #Descomentar si se desea comenzar el script desde este punto (no olvidar cargar paquetes)
datos_03
## # A tibble: 255 × 17
## Hexágono Código Coordenadas Datum Altura Incert…¹ Especie Colec…²
## <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 12 H12M1 18.4675906 -69.9139862 WGS84 53 msnm 5 metros Pheido… D. Guz…
## 2 13 H13M2 18.4672803 -69.9192371 WGS84 54 msnm 4 metros Dorymy… D. Guz…
## 3 13 H13M3 18.4672054 -69.9192181 WGS84 54 msnm 5 metros Paratr… D. Guz…
## 4 13 H13M4 18.4669508 -69.9193136 WGS84 54 msnm 7 metros Pheido… D. Guz…
## 5 13 H13M5 18.4671236 -69.9191236 WGS84 54 msnm 3 metros Soleno… D. Guz…
## 6 13 H13M6 18.4671876 -69.9193009 WGS84 53 msnm 4 metros Pheido… D. Guz…
## 7 12 H12M7 18.4675406 -69.9194121 WGS84 53 msnm 7 metros Pheido… D. Guz…
## 8 12 H12M8 18.4675806 -69.9193971 WGS84 53 msnm 5 metros Pheido… D. Guz…
## 9 12 H12M9 18.4675541 -69.9193724 WGS84 53 msnm 7 metros Pheido… D. Guz…
## 10 12 H12M10 18.4678106 -69.9148187 WGS84 53 msnm 7 metros Soleno… D. Guz…
## # … with 245 more rows, 9 more variables: Fecha <dttm>,
## # `No. de individos` <dbl>, Observaciones <chr>, País <chr>, Provincia <chr>,
## # Municipio <chr>, `Área protegida` <chr>, lat <dbl>, lon <dbl>, and
## # abbreviated variable names ¹Incertidumbre, ²`Colector/os`
generar_mc <- function(patron = '^H.*', presencia_ausencia = T){
datos_03 %>%
filter(grepl(patron, Código)) %>%
select(Hexágono, Especie) %>%
mutate(n = 1) %>%
pivot_wider(names_from = Especie, values_from = n,
values_fn = ~ sum(.x, na.rm = TRUE),
values_fill = 0) %>%
arrange(Hexágono) %>%
column_to_rownames('Hexágono') %>%
{if(presencia_ausencia) replace(., . > 1, 1) else .}
}
mc_marzo <- generar_mc(patron = '^H.*')
# png('grafico_de_mosaicos_marzo.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(
add_row(mc_marzo, .before = 1) %>% replace(is.na(.), 0) %>% remove_rownames()
)
# dev.off()
mc_mayo <- generar_mc(patron = '^2H.*')
# png('grafico_de_mosaicos_mayo.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(mc_mayo)
# dev.off()
mc_pooled <- generar_mc(patron = '.*')
# png('grafico_de_mosaicos_pooled.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(mc_pooled)
# dev.off()
mc_marzo_nidos <- generar_mc(patron = '^H.*', presencia_ausencia = F)
# png('grafico_de_mosaicos_nidos_marzo.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(
add_row(mc_marzo_nidos, .before = 1) %>% replace(is.na(.), 0) %>% remove_rownames()
)
# dev.off()
mc_mayo_nidos <- generar_mc(patron = '^2H.*', presencia_ausencia = F)
# png('grafico_de_mosaicos_nidos_mayo.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(mc_mayo_nidos)
# dev.off()
# Nidos combinados
mc_nidos <- generar_mc(patron = '.*', presencia_ausencia = F)
# png('grafico_de_mosaicos_nidos.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(mc_nidos)
# dev.off()
# Reescalado
mc_puente <- scales::rescale(as.matrix(mc_nidos), to= c(1, 5))
mc_nidos_reescalado <- round(ifelse(mc_puente==1, 0, mc_puente), 0)
mc_todas <- list(mc_marzo, mc_mayo, mc_pooled)
names(mc_todas) <- c('Marzo', 'Mayo', 'Combinada')
#' ### Lista de especies
map(mc_todas,
~ sort(colnames(.x)) %>%
as_tibble() %>%
mutate(fila = 1:nrow(.)) %>%
select(2, 1) %>%
kable(col.names = c('#', 'Especie'), table.attr='class="myTable"'))
$Marzo
| # | Especie |
|---|---|
| 1 | Cyphomyrmex rimosus |
| 2 | Dorymyrmex antillanus |
| 3 | Mycetomoellerius jamaicensis |
| 4 | Odontomachus ruginodis |
| 5 | Paratrechina longicornis |
| 6 | Pheidole jelskii |
| 7 | Pheidole moerens |
| 8 | Pheidole subarmata |
| 9 | Solenopsis geminata |
| 10 | Tetramorium bicarinatum |
| 11 | Tetramorium lanuginosum |
| 12 | Wasmannia auropunctata |
| # | Especie |
|---|---|
| 1 | Acropyga dubitata |
| 2 | Brachymyrmex heeri |
| 3 | Cyphomyrmex rimosus |
| 4 | Dorymyrmex antillanus |
| 5 | Monomorium floricola |
| 6 | Mycetomoellerius jamaicensis |
| 7 | Odontomachus ruginodis |
| 8 | Paratrechina longicornis |
| 9 | Pheidole jelskii |
| 10 | Pheidole subarmata |
| 11 | Pogonomyrmex schmitti |
| 12 | Solenopsis geminata |
| 13 | Tetramorium lanuginosum |
| 14 | Wasmannia auropunctata |
| # | Especie |
|---|---|
| 1 | Acropyga dubitata |
| 2 | Brachymyrmex heeri |
| 3 | Cyphomyrmex rimosus |
| 4 | Dorymyrmex antillanus |
| 5 | Monomorium floricola |
| 6 | Mycetomoellerius jamaicensis |
| 7 | Odontomachus ruginodis |
| 8 | Paratrechina longicornis |
| 9 | Pheidole jelskii |
| 10 | Pheidole moerens |
| 11 | Pheidole subarmata |
| 12 | Pogonomyrmex schmitti |
| 13 | Solenopsis geminata |
| 14 | Tetramorium bicarinatum |
| 15 | Tetramorium lanuginosum |
| 16 | Wasmannia auropunctata |
cobertura <- st_read('coberturas.gpkg')
## Reading layer `coberturas' from data source
## `C:\Users\User\Desktop\Universidad\Tesis\Datos de los muestreos\coberturas.gpkg'
## using driver `GPKG'
## Simple feature collection with 82 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 402699.2 ymin: 2041878 xmax: 403116 ymax: 2042246
## Projected CRS: WGS 84 / UTM zone 19N
hex_inters <- st_as_sf(st_intersection(st_union(cobertura), hex)) %>% mutate(ENLACE = hex$ENLACE)
ma_sf <- zstatsf(zones = hex_inters, values = cobertura, grpx = "ENLACE", grpy = "tipo")
ma_sf <- ma_sf %>% mutate(across(2:5, units::drop_units))
ma_sf %>% plot()
ma <- ma_sf %>% st_drop_geometry() %>% column_to_rownames('ENLACE')
ma
## construido, mobiliario (bordes edificios, acerado, bancos, postes...)
## 1 0.0000000
## 2 18.4869756
## 3 0.0000000
## 4 3.4888882
## 5 0.0000000
## 6 2.4769313
## 7 0.6085964
## 8 27.2501334
## 9 7.5820449
## 10 0.1769173
## 11 5.1515836
## 12 38.0738217
## 13 51.8993681
## 14 0.0000000
## 15 1.1269871
## 16 3.7630314
## 17 69.1243002
## 18 16.2492695
## 19 25.6929698
## 20 15.8724897
## 21 21.0199767
## 22 43.0273268
## 23 11.3836316
## 24 44.3322007
## 25 23.9599716
## 26 1.7747530
## 27 18.4280883
## 28 0.3534807
## dosel edificación erguida suelo, herbáceas, no edificado ni cubierto
## 1 89.86997 0.000000 10.13003022
## 2 80.17277 0.000000 1.34025095
## 3 99.88796 0.000000 0.11204383
## 4 61.57379 24.838187 10.09913880
## 5 77.72845 22.271546 0.00000000
## 6 95.53464 0.000000 1.98843129
## 7 97.26926 0.000000 2.12214164
## 8 67.74032 0.000000 5.00955070
## 9 54.37792 37.902586 0.13744975
## 10 88.68445 0.000000 11.13863183
## 11 80.03134 1.765534 13.05153872
## 12 52.72902 0.000000 9.19716130
## 13 46.72021 1.068254 0.31216408
## 14 96.50561 0.000000 3.49438827
## 15 97.93200 0.000000 0.94101640
## 16 96.23697 0.000000 0.00000000
## 17 26.98363 0.000000 3.89206487
## 18 63.99899 19.294636 0.45710407
## 19 74.29148 0.000000 0.01555076
## 20 84.12751 0.000000 0.00000000
## 21 67.26251 0.000000 11.71751267
## 22 15.30474 28.156893 13.51103755
## 23 41.44998 46.221870 0.94452067
## 24 53.38524 0.000000 2.28256075
## 25 76.04003 0.000000 0.00000000
## 26 65.53121 24.489490 8.20454626
## 27 78.54138 0.000000 3.03053153
## 28 93.44066 0.000000 6.20586141
rowSums(ma) #Comprobación de suma 100 por filas
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
## 21 22 23 24 25 26 27 28
## 100 100 100 100 100 100 100 100
#' ### Número de sitios, tanto en matriz de comunidad como en ambiental
#' Verifica que coinciden
nrow(mc_todas$Combinada) #En la matriz de comunidad
## [1] 28
nrow(ma) #En la matriz ambiental
## [1] 28
#' ### Riqueza numérica de especies (usando matriz de comunidad) por quadrat
#' Nota: cargar paquete vegan arriba, en el área de paquetes
specnumber(mc_todas$Combinada)
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
## 2 3 2 4 4 7 9 3 4 2 6 5 6 2 6 3 1 3 4 4 6 3 2 4 5 2
## 27 28
## 4 4
sort(specnumber(mc_todas$Combinada)) # Ordenados ascendentemente
## 17 1 3 10 14 23 26 2 8 16 18 22 4 5 9 19 20 24 27 28 12 25 11 13 15 21
## 1 2 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 4 4 4 5 5 6 6 6 6
## 6 7
## 7 9
summary(specnumber(mc_todas$Combinada)) # Resumen estadístico
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.000 2.750 4.000 3.929 5.000 9.000
#' ### Riqueza numérica de toda la "comunidad"
specnumber(colSums(mc_todas$Combinada))
## [1] 16
En esta seccion se realiza una clasificacion no supervisada de habitats en funcion de su composicion de especies de hormigas, y se explora cuales atributos (de los 4 generados a partir del mapa de cobertura) caracterizan a dichos habitats.
mc_orig <- mc_todas$Combinada
data.frame(`Número de hexágonos` = sort(colSums(mc_orig), decreasing = T), check.names = F) %>%
kable(booktabs=T) %>%
kable_styling(latex_options = c("HOLD_position", "scale_down")) %>%
gsub(' NA ', '', .) # Número de hexágonos en los que está presente cada especie
| Número de hexágonos | |
|---|---|
| Pheidole subarmata | 25 |
| Solenopsis geminata | 20 |
| Dorymyrmex antillanus | 19 |
| Paratrechina longicornis | 8 |
| Pheidole jelskii | 8 |
| Wasmannia auropunctata | 6 |
| Cyphomyrmex rimosus | 6 |
| Odontomachus ruginodis | 5 |
| Tetramorium lanuginosum | 3 |
| Pogonomyrmex schmitti | 3 |
| Mycetomoellerius jamaicensis | 2 |
| Pheidole moerens | 1 |
| Tetramorium bicarinatum | 1 |
| Brachymyrmex heeri | 1 |
| Monomorium floricola | 1 |
| Acropyga dubitata | 1 |
# Usa el vector anterior para determinar un umbral o rango de registros para filtrar tu matriz
# ¿En cuántos hexágonos está cada especie? Filtra tus datos usando tu propio criterio.
# Especies que aparecen en pocos hexágonos se consideran "raras". Por ejemplo, si una especie sólo
# aparece en un hexágono en todo el país, es un "singleton", si en dos, "doubleton", y así.
# Estas especies podrían contribuir a generar "ruido" en análisis posteriores, se recomienda excluirlas.
# Elige un valor mínimo (representado por único número entero) o por un rango de enteros (e.g. de 10 a 20),
# para seleccionar las especies que estén mejor representadas de acuerdo a tu criterio.
# Por ejemplo, si usas el valor m, el script considerará a este valor como "el número mínimo de hexágonos
# en los que está representada una especie, y creará una matriz de comunidad de especies seleccionadas
# que están presentes en m hexágonos o más. Si eliges un rango, por ejemplo [m,n], el script generará
# una matriz de comunidad que representadas un mínimo de m hexágonos y un máximo de n hexágonos.
# (ambos extremos inclusive).
en_cuantos_hex <- 2 #!!!ESTE UMBRAL EXCLUYE A ESPECIES QUE APARECEN EN SOLO 1 HEXAGONO, Y CONSERVA SOLO AQUELLAS QUE APARECEN EN 2 O MAS HEXAGONOS
# Explicación: "en_cuantos_hex <- X", donde X es el número de hexágonos mínimo donde cada especie
# debe estar presente. IMPORTANTE: elige TU PROPIO umbral.
{if(length(en_cuantos_hex)==1) selector <- en_cuantos_hex:max(colSums(mc_orig)) else
if(length(en_cuantos_hex)==2)
selector <- min(en_cuantos_hex):max(en_cuantos_hex) else
stop('Debes indicar uno o dos valores numéricos')}
selector
## [1] 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
mc_orig_seleccionadas <- mc_orig[, colSums(mc_orig) %in% selector]
# Mínimo número de especies por hexágono
data.frame(`Número de especies por hexágono` = sort(rowSums(mc_orig), decreasing = T), check.names = F) %>%
kable(booktabs=T) %>%
kable_styling(latex_options = c("HOLD_position", "scale_down")) %>%
gsub(' NA ', '', .) # Número de hexágonos en los que está presente cada especie
| Número de especies por hexágono | |
|---|---|
| 7 | 9 |
| 6 | 7 |
| 11 | 6 |
| 13 | 6 |
| 15 | 6 |
| 21 | 6 |
| 12 | 5 |
| 25 | 5 |
| 4 | 4 |
| 5 | 4 |
| 9 | 4 |
| 19 | 4 |
| 20 | 4 |
| 24 | 4 |
| 27 | 4 |
| 28 | 4 |
| 2 | 3 |
| 8 | 3 |
| 16 | 3 |
| 18 | 3 |
| 22 | 3 |
| 1 | 2 |
| 3 | 2 |
| 10 | 2 |
| 14 | 2 |
| 23 | 2 |
| 26 | 2 |
| 17 | 1 |
min_especies_por_hex <- 2 #!!!ESTA LINEA SACA LOS HEXAGONOS QUE TIENEN SOLO 1 ESPECIE, Y CONSERVA SOLO AQUELLOS QUE TIENEN 2 O MAS ESPECIES PRESENTES (EN EL CASO DE LA MATRIZ DE ESPECIES COMBINADAS, EXCLUYE EL HEXAGONO 17)
# Explicación: "min_especies_por_hex <- Y", donde Y es el número mínimo (inclusive) de especies
# que debe existir en cada hexágono. Por debajo de dicho valor, el hexágono es excluido.
mi_fam <- mc_orig_seleccionadas[rowSums(mc_orig_seleccionadas)>=min_especies_por_hex, ]
nrow(mi_fam)
## [1] 27
# mi_fam <- mc_orig_seleccionadas[!rowSums(mc_orig_seleccionadas)==0, ] #Elimina filas sin registros
# rowSums(mi_fam) #Riqueza por hexágonos con especies seleccionadas. Comentado por extenso
all(rowSums(mi_fam)>0) #Debe ser TRUE: todos los hexágonos tienen al menos 1 registro
## [1] TRUE
ncol(mi_fam) #Riqueza de especies
## [1] 11
# IMPORTANTE, quitar Formicidae. ¡¡Muchos nombres que habría que resolver si se quisiera publicar!!
# Usar nombres cortos o abreviados para las especies
# nombres_largos <- colnames(mi_fam)
# (colnames(mi_fam) <- make.cepnames(word(colnames(mi_fam), 1, 2)))
# (df_equivalencias <- data.frame(
# nombre_original = nombres_largos,
# abreviado = colnames(mi_fam)))
En esta subseccion se realiza propiamente el analisis de agrupamiento segun cuatro metodos distintos (single, completate, upgma, Ward); se elige, al final, el metodo Ward. Se contrastan los metodos con remuestreo multiescalas por bootstrap.
mi_fam_t <- decostand(mi_fam, 'hellinger') #Hellinger
env <- ma[match(rownames(mi_fam), rownames(ma)), ]
all(rownames(mi_fam) == rownames(env)) #Si es TRUE, sigue adelante
## [1] TRUE
library(broom)
library(cluster)
library(gclus)
library(pvclust)
mi_fam_d <- vegdist(mi_fam_t, "euc")
lista_cl <- list(
cl_single = hclust(mi_fam_d, method = 'single'),
cl_complete = hclust(mi_fam_d, method = 'complete'),
cl_upgma = hclust(mi_fam_d, method = 'average'),
cl_ward = hclust(mi_fam_d, method = 'ward.D2')
)
par(mfrow = c(2,2))
invisible(map(names(lista_cl), function(x) plot(lista_cl[[x]], main = x, hang = -1, cex = 1)))
par(mfrow = c(1,1))
cutree(lista_cl$cl_ward, 2)
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 18 19 20 21 22 23 24 25 26 27
## 1 1 2 1 1 2 2 2 2 1 2 1 1 1 2 1 1 1 1 1 1 1 1 2 2 1
## 28
## 2
if(interactive()) dev.new()
cl_pvclust_ward <-
pvclust(t(mi_fam_t),
method.hclust = "ward.D2",
method.dist = "euc",
iseed = 999, # Resultado reproducible
parallel = TRUE)
## Creating a temporary cluster...done:
## socket cluster with 3 nodes on host 'localhost'
## Multiscale bootstrap... Done.
# Añadir los valores de p
plot(cl_pvclust_ward, hang = -1)
# Añadir rectángulos a los grupos significativos
lines(cl_pvclust_ward)
pvrect(cl_pvclust_ward, alpha = 0.91, border = 4)
w_dend_reord <- reorder.hclust(lista_cl$cl_ward, mi_fam_d)
grupos_ward_k2 <- as.factor(cutree(lista_cl$cl_ward, k = 2))
table(grupos_ward_k2) #¿Cuántos hexágonos hay en cada grupo?
## grupos_ward_k2
## 1 2
## 17 10
plot(w_dend_reord, hang = -1)
rect.hclust(tree = w_dend_reord, k = 2)
grupos_ward_k3 <- as.factor(cutree(lista_cl$cl_ward, k = 3))
table(grupos_ward_k3) #¿Cuántos hexágonos hay en cada grupo?
## grupos_ward_k3
## 1 2 3
## 6 10 11
plot(w_dend_reord, hang = -1)
rect.hclust(tree = w_dend_reord, k = 3)
# Guardaré estos vectores en archivos para reutilizarlos en *scripts* posteriores:
# saveRDS(grupos_ward_k2, 'grupos_ward_k2.RDS')
# saveRDS(grupos_ward_k3, 'grupos_ward_k3.RDS') #No hubo patrones significativos de asociacion
ma_para_combinadas <- ma[match(rownames(mi_fam), rownames(ma)), ]
(m_amb_ward_k2 <- ma_para_combinadas %>%
rownames_to_column('hex_id') %>%
mutate(grupos_ward_k2) %>%
pivot_longer(-c(grupos_ward_k2, hex_id), names_to = "variable", values_to = "valor"))
## # A tibble: 108 × 4
## hex_id grupos_ward_k2 variable valor
## <chr> <fct> <chr> <dbl>
## 1 1 1 construido, mobiliario (bordes edificios, acerad… 0
## 2 1 1 dosel 89.9
## 3 1 1 edificación erguida 0
## 4 1 1 suelo, herbáceas, no edificado ni cubierto 10.1
## 5 2 1 construido, mobiliario (bordes edificios, acerad… 18.5
## 6 2 1 dosel 80.2
## 7 2 1 edificación erguida 0
## 8 2 1 suelo, herbáceas, no edificado ni cubierto 1.34
## 9 3 2 construido, mobiliario (bordes edificios, acerad… 0
## 10 3 2 dosel 99.9
## # … with 98 more rows
m_amb_ward_k2_tw <- m_amb_ward_k2 %>%
group_by(variable) %>%
summarise(
p_valor_t = tryCatch(t.test(valor ~ grupos_ward_k2)$p.value, error = function(e) NA),
p_valor_w = tryCatch(wilcox.test(valor ~ grupos_ward_k2)$p.value, error = function(e) NA)
) %>%
arrange(p_valor_t)
m_amb_ward_k2_tw %>%
kable(booktabs=T) %>%
kable_styling(latex_options = c("HOLD_position", "scale_down")) %>%
gsub(' NA ', '', .)
| variable | p_valor_t | p_valor_w |
|---|---|---|
| construido, mobiliario (bordes edificios, acerado, bancos, postes…) | 0.0405413 | 0.1910607 |
| dosel | 0.0696707 | 0.0743239 |
| suelo, herbáceas, no edificado ni cubierto | 0.6704152 | 0.8405846 |
| edificación erguida | 0.7283481 | 0.8108505 |
m_amb_ward_k2 %>%
group_by(variable) %>%
ggplot() + aes(x = grupos_ward_k2, y = valor, fill = grupos_ward_k2) +
geom_boxplot(lwd = 0.2) +
scale_fill_brewer(palette = 'Set1') +
theme_bw(base_size=12) +
theme(legend.position="none") +
facet_wrap(~ variable, scales = 'free_y', ncol = 2,
labeller = label_wrap_gen(width = 35,multi_line = T))
# CON WARD K3 NO APARECEN EFECTOS SIGNIFICATIVOS. MEJOR USAR SOLO AGRUPAMIENTO WARD K2
En resumen, las variables de porcentaje de “construido…” y “dosel…” son las unicas que presentaron efectos significativos entre grupos, por lo que consideramos que estas son las que, con mayor probabilidad, se asocian con la composicion de especies diferenciada entre grupos. El grupo 1 se caracteriza por tener areas construidas predominantemente, y el 2 por tener predominio de dosel.
phi_ward_k2 <- multipatt(
mi_fam,
grupos_ward_k2,
func = "r.g",
max.order = 1,
control = how(nperm = 999))
summary(phi_ward_k2)
##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.05
##
## Total number of species: 11
## Selected number of species: 3
## Number of species associated to 1 group: 3
##
## List of species associated to each combination:
##
## Group 1 #sps. 1
## stat p.value
## Dorymyrmex antillanus 0.749 0.001 ***
##
## Group 2 #sps. 2
## stat p.value
## Wasmannia auropunctata 0.655 0.002 **
## Cyphomyrmex rimosus 0.492 0.016 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ma_sf_grupos_ward_k2 <- ma_sf %>%
inner_join(m_amb_ward_k2 %>%
select(ENLACE=hex_id, grupos_ward_k2) %>%
distinct() %>%
mutate(ENLACE=as.integer(ENLACE)))
st_geometry(ma_sf_grupos_ward_k2) <- 'geometry'
mapa_ma_sf_grupos_ward_k2 <- ma_sf_grupos_ward_k2 %>%
ggplot + aes(fill = grupos_ward_k2) +
geom_sf()
if(interactive()) dev.new()
mapa_ma_sf_grupos_ward_k2
mapa_ma_sf <- ma_sf_grupos_ward_k2 %>%
select(-grupos_ward_k2) %>%
pivot_longer(-c(geometry, ENLACE)) %>%
group_by(name) %>%
mutate(value=scale(value)) %>%
ggplot + aes(fill = value) +
scale_fill_viridis_b() +
geom_sf() +
facet_wrap(~name, ncol = 2)
if(interactive()) dev.new()
mapa_ma_sf
Analisis de ordenacion. No aporta patrones de interes, puesto que la inercia restringida es muy pequena (es habitual cuando se trabaja con datos de presencia/ausencia). Analisis descartado en este caso.
# mi_fam_t_rda <- rda(mi_fam_t ~ .,
# ma[-17,] %>% select(all_of(matches('construido|dosel'))),
# scale = T)
# summary(mi_fam_t_rda)
# RsquareAdj(mi_fam_t_rda)$adj.r.squared
# vif.cca(mi_fam_t_rda)
# escalado <- 1
# plot(mi_fam_t_rda,
# scaling = escalado,
# display = c("sp", "lc", "cn"),
# main = paste("Triplot de RDA especies ~ var. GSL + ESA, escalamiento", escalado)
# )
# mi_fam_t_rda_sc1 <- scores(mi_fam_t_rda,
# choices = 1:2,
# scaling = escalado,
# display = "sp"
# )
# # text(mi_fam_t_rda, "species", col="red", cex=0.8, scaling=escalado)
# arrows(0, 0,
# mi_fam_t_rda_sc1[, 1] * 0.9,
# mi_fam_t_rda_sc1[, 2] * 0.9,
# length = 0,
# lty = 1,
# col = "red"
# )
Pool de especies, estimadores tradicionales de la riqueza (Chao, Jack1, Jack2, etc.). Esta matriz tiene muchas especies raras, por lo que los estimadores penalizan la riqueza capturada. Se considerara usar una matriz ponderada por el numero de nidos.
# mc_orig <- mc_todas$Combinada
specpool(mc_orig)
## Species chao chao.se jack1 jack1.se jack2 boot boot.se n
## All 16 28.05357 16.54004 20.82143 2.564723 24.57011 18.0219 1.375846 28
specpool(mc_orig)[[1]]/specpool(mc_orig)[-c(3,5,8)]*100
## Species chao jack1 jack2 boot n
## All 100 57.03374 76.84391 65.11978 88.78088 57.14286
# Nidos
estimateR(colSums(mc_nidos))
## S.obs S.chao1 se.chao1 S.ACE se.ACE
## 16.000000 26.000000 10.258674 20.465909 2.201645
estimateR(colSums(mc_nidos))[[1]]/estimateR(colSums(mc_nidos))[-c(3,5,8)]*100
## S.obs S.chao1 S.ACE
## 100.00000 61.53846 78.17879
Estimadores
Matriz original combinada: el error “Error in if (var_iChao2 > 0) { : missing value where TRUE/FALSE needed” impide ejecutar la estimacion para la matriz combinada.
# Incidencia
ChaoSpecies(data.frame(V1 = c(nrow(mc_orig), as.integer(colSums(mc_orig)))),
datatype = 'incidence_freq', k=1, conf=0.95)
# Nmero de nidos
ChaoSpecies(data.frame(V1 = c(nrow(mc_nidos), as.integer(colSums(mc_nidos)))),
datatype = 'abundance', k=1, conf=0.95)
## Warning: In this case, it can't estimate the variance of Homogeneous estimation
##
## (1) BASIC DATA INFORMATION:
##
## Variable Value
## Sample size n 283
## Number of observed species D 17
## Coverage estimate for entire dataset C 0.982
## CV for entire dataset CV 1.538
## Cut-off point k 1
##
## Variable Value
## Number of observed individuals for rare group n_rare 5
## Number of observed species for rare group D_rare 5
## Estimate of the sample coverage for rare group C_rare 0
## Estimate of CV for rare group in ACE CV_rare 0
## Estimate of CV1 for rare group in ACE-1 CV1_rare 0
## Number of observed individuals for abundant group n_abun 278
## Number of observed species for abundant group D_abun 12
##
## NULL
##
##
## (2) SPECIES RICHNESS ESTIMATORS TABLE:
##
## Estimate s.e. 95%Lower 95%Upper
## Homogeneous Model Inf NA NA NA
## Homogeneous (MLE) 17.000 1.123 17.000 21.705
## Chao1 (Chao, 1984) 26.965 10.235 18.888 69.583
## Chao1-bc 26.965 10.235 18.888 69.583
## iChao1 (Chiu et al. 2014) 29.465 10.607 19.934 69.955
## ACE (Chao & Lee, 1992) 26.965 10.235 18.888 69.583
## ACE-1 (Chao & Lee, 1992) 26.965 10.235 18.888 69.583
## 1st order jackknife 21.982 3.154 18.597 32.544
## 2nd order jackknife 26.947 5.453 20.642 44.166
##
##
## (3) DESCRIPTION OF ESTIMATORS/MODELS:
##
## Homogeneous Model: This model assumes that all species have the same incidence or detection probabilities. See Eq. (3.2) of Lee and Chao (1994) or Eq. (12a) in Chao and Chiu (2016b).
##
## Chao2 (Chao, 1987): This approach uses the frequencies of uniques and duplicates to estimate the number of undetected species; see Chao (1987) or Eq. (11a) in Chao and Chiu (2016b).
##
## Chao2-bc: A bias-corrected form for the Chao2 estimator; see Chao (2005).
##
## iChao2: An improved Chao2 estimator; see Chiu et al. (2014).
##
## ICE (Incidence-based Coverage Estimator): A non-parametric estimator originally proposed by Lee and Chao (1994) in the context of capture-recapture data analysis. The observed species are separated as frequent and infrequent species groups; only data in the infrequent group are used to estimate the number of undetected species. The estimated CV for species in the infrequent group characterizes the degree of heterogeneity among species incidence probabilities. See Eq. (12b) of Chao and Chiu (2016b), which is an improved version of Eq. (3.18) in Lee and Chao (1994). This model is also called Model(h) in capture-recapture literature where h denotes "heterogeneity".
##
## ICE-1: A modified ICE for highly-heterogeneous cases.
##
## 1st order jackknife: It uses the frequency of uniques to estimate the number of undetected species; see Burnham and Overton (1978).
##
## 2nd order jackknife: It uses the frequencies of uniques and duplicates to estimate the number of undetected species; see Burnham and Overton (1978).
##
## 95% Confidence interval: A log-transformation is used for all estimators so that the lower bound of the resulting interval is at least the number of observed species. See Chao (1987).
Para fines didacticos, representamos la curva de acumulacion de especies de la matriz combinada, donde se observa una pronunciada pendiente y una escasa capacidad de alcanzarse la asintota. Decimos “para fines didacticos”, porque complementa el analisis posterior por grupos segun Ward.
mc_orig_lista <- sapply(rownames(mc_orig), function(x) as.numeric(mc_orig[x,]), simplify = F)
nasin_raref_pooled <- iNEXT::iNEXT(
x = colSums(mc_orig), #if(is.list(mc_orig)) mc_orig_lista else mc_orig,
q=0,
knots = 400,
datatype = 'incidence_freq')
acumulacion_especies_pooled <- iNEXT::ggiNEXT(nasin_raref_pooled, type=1) +
theme_bw() +
theme(
text = element_text(size = 20),
panel.background = element_rect(fill = 'white', colour = 'black'),
panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
) +
ylab('Riqueza de especies') +
xlab('Número de sitios') +
scale_y_continuous(breaks = seq(0,80, length.out = 9)) +
scale_color_manual(values = brewer.pal(8, 'Set2')) +
scale_fill_manual(values = brewer.pal(8, 'Set2'))
if(interactive()) dev.new()
acumulacion_especies_pooled
Numero de nidos
mc_nidos_lista <- sapply(rownames(mc_nidos), function(x) as.numeric(mc_nidos[x,]), simplify = F)
nasin_raref_nidos <- iNEXT::iNEXT(
x = colSums(mc_nidos), #if(is.list(mc_nidos)) mc_nidos_lista else mc_nidos,
q=0,
knots = 400,
datatype = 'abundance')
acumulacion_especies_nidos <- iNEXT::ggiNEXT(nasin_raref_nidos, type=1) +
theme_bw() +
theme(
text = element_text(size = 20),
panel.background = element_rect(fill = 'white', colour = 'black'),
panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
) +
ylab('Riqueza de especies') +
xlab('Número de sitios') +
scale_y_continuous(breaks = seq(0,80, length.out = 9)) +
scale_color_manual(values = brewer.pal(8, 'Set2')) +
scale_fill_manual(values = brewer.pal(8, 'Set2'))
if(interactive()) dev.new()
acumulacion_especies_nidos
Numero de nidos reescalado
mc_nidos_reescalado_lista <- sapply(rownames(mc_nidos_reescalado), function(x) as.numeric(mc_nidos_reescalado[x,]), simplify = F)
nasin_raref_nidos_reescalado <- iNEXT::iNEXT(
x = colSums(mc_nidos_reescalado), #if(is.list(mc_nidos)) mc_nidos_lista else mc_nidos,
q=0,
knots = 400,
datatype = 'abundance')
acumulacion_especies_nidos_reescalado <- iNEXT::ggiNEXT(nasin_raref_nidos_reescalado, type=1) +
theme_bw() +
theme(
text = element_text(size = 20),
panel.background = element_rect(fill = 'white', colour = 'black'),
panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
) +
ylab('Riqueza de especies') +
xlab('Número de sitios') +
scale_y_continuous(breaks = seq(0,80, length.out = 9)) +
scale_color_manual(values = brewer.pal(8, 'Set2')) +
scale_fill_manual(values = brewer.pal(8, 'Set2'))
if(interactive()) dev.new()
acumulacion_especies_nidos_reescalado
Se refiere a la matriz filtrada donde se eliminaron especies poco representadas
specpool(mi_fam)
## Species chao chao.se jack1 jack1.se jack2 boot boot.se n
## All 11 11 0 11 0 10.10969 11.21473 0.4439909 27
ChaoSpecies(data.frame(V1 = c(nrow(mi_fam), as.numeric(colSums(mi_fam)))),
datatype = 'incidence_freq', k=2, conf=0.95)
## Warning: In this case, it can't estimate the variance of Homogeneous estimation
##
## Warning: In this case, it can't estimate the variance of Model(h) estimation
##
## Warning: In this case, it can't estimate the variance of Model(h)-1 estimation
##
## Warning: In this case, it can't estimate the variance of 1st-order-jackknife estimation
##
## Warning: In this case, it can't estimate the variance of 2nd-order-jackknife estimation
##
## Warning: In this case, it can't estimate the variance of Model(h) estimation
##
## Warning: In this case, it can't estimate the variance of Model(h)-1 estimation
##
## (1) BASIC DATA INFORMATION:
##
## Variable Value
## Number of observed species D 11
## Number of sampling units T 27
## Total number of incidences U 104
## Coverage estimate for entire dataset C 1
## CV for entire dataset CV 0.759
##
## Variable Value
## Cut-off point k 2
## Total number of incidences in infrequent group U_infreq 2
## Number of observed species for infrequent group D_infreq 1
## Estimated sample coverage for infrequent group C_infreq 1
## Estimated CV for infrequent group in ICE CV_infreq 0.196
## Estimated CV1 for infrequent group in ICE-1 CV1_infreq 0.196
## Number of observed species for frequent group D_freq 10
##
## Q1 Q2
## Incidence freq. counts 0 1
##
##
## (2) SPECIES RICHNESS ESTIMATORS TABLE:
##
## Estimate s.e. 95%Lower 95%Upper
## Homogeneous Model 11.00 0.457 11 12.169
## Chao2 (Chao, 1987) 11.00 0.457 11 12.169
## Chao2-bc 11.00 0.457 11 12.169
## iChao2 (Chiu et al. 2014) 11.00 0.457 11 12.169
## ICE (Lee & Chao, 1994) 11.00 0.457 11 12.169
## ICE-1 (Lee & Chao, 1994) 11.00 0.457 11 12.169
## 1st order jackknife 11.00 0.457 11 12.169
## 2nd order jackknife 10.11 NA NA NA
##
##
## (3) DESCRIPTION OF ESTIMATORS/MODELS:
##
## Homogeneous Model: This model assumes that all species have the same incidence or detection probabilities. See Eq. (3.2) of Lee and Chao (1994) or Eq. (12a) in Chao and Chiu (2016b).
##
## Chao2 (Chao, 1987): This approach uses the frequencies of uniques and duplicates to estimate the number of undetected species; see Chao (1987) or Eq. (11a) in Chao and Chiu (2016b).
##
## Chao2-bc: A bias-corrected form for the Chao2 estimator; see Chao (2005).
##
## iChao2: An improved Chao2 estimator; see Chiu et al. (2014).
##
## ICE (Incidence-based Coverage Estimator): A non-parametric estimator originally proposed by Lee and Chao (1994) in the context of capture-recapture data analysis. The observed species are separated as frequent and infrequent species groups; only data in the infrequent group are used to estimate the number of undetected species. The estimated CV for species in the infrequent group characterizes the degree of heterogeneity among species incidence probabilities. See Eq. (12b) of Chao and Chiu (2016b), which is an improved version of Eq. (3.18) in Lee and Chao (1994). This model is also called Model(h) in capture-recapture literature where h denotes "heterogeneity".
##
## ICE-1: A modified ICE for highly-heterogeneous cases.
##
## 1st order jackknife: It uses the frequency of uniques to estimate the number of undetected species; see Burnham and Overton (1978).
##
## 2nd order jackknife: It uses the frequencies of uniques and duplicates to estimate the number of undetected species; see Burnham and Overton (1978).
##
## 95% Confidence interval: A log-transformation is used for all estimators so that the lower bound of the resulting interval is at least the number of observed species. See Chao (1987).
Curva de acumulacion de especies. Primero generamos la matriz combinada (pooled) por grupos del agrupamiento Ward K2.
mi_fam_ward <- mi_fam %>%
mutate(g = grupos_ward_k2) %>%
group_by(g) %>%
summarise_all(sum) %>%
select(-g) %>%
mutate(N = nrow(mi_fam)) %>%
relocate(N, .before = 1) %>%
data.frame
mi_fam_ward
## N Pheidole.subarmata Dorymyrmex.antillanus Paratrechina.longicornis
## 1 27 16 16 3
## 2 27 9 2 5
## Solenopsis.geminata Pheidole.jelskii Mycetomoellerius.jamaicensis
## 1 11 6 1
## 2 9 2 1
## Tetramorium.lanuginosum Odontomachus.ruginodis Wasmannia.auropunctata
## 1 1 2 0
## 2 2 3 6
## Cyphomyrmex.rimosus Pogonomyrmex.schmitti
## 1 1 1
## 2 5 2
Posteriormente, generamos la curva de acumlacion.
nasin_raref <- iNEXT::iNEXT(
x = t(mi_fam_ward),
q=0,
knots = 400,
datatype = 'incidence_freq')
acumulacion_especies <- iNEXT::ggiNEXT(nasin_raref, type=1) +
theme_bw() +
theme(
text = element_text(size = 20),
panel.background = element_rect(fill = 'white', colour = 'black'),
panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
) +
ylab('Riqueza de especies') +
xlab('Número de sitios') +
scale_y_continuous(breaks = seq(0,80, length.out = 9)) +
scale_color_manual(values = brewer.pal(8, 'Set2')) +
scale_fill_manual(values = brewer.pal(8, 'Set2'))
if(interactive()) dev.new()
acumulacion_especies
#' ### Número de sitios, tanto en matriz de comunidad como en ambiental
#' Verifica que coinciden
nrow(mc_todas$Marzo) #En la matriz de comunidad
## [1] 27
nrow(ma) #En la matriz ambiental
## [1] 28
#' ### Riqueza numérica de especies (usando matriz de comunidad) por quadrat
#' Nota: cargar paquete vegan arriba, en el área de paquetes
specnumber(mc_todas$Marzo)
## 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
## 2 1 3 3 2 2 2 2 1 2 5 5 1 2 2 1 3 3 3 4 2 1 2 4 2 2
## 28
## 1
sort(specnumber(mc_todas$Marzo)) # Ordenados ascendentemente
## 3 10 14 17 23 28 2 6 7 8 9 11 15 16 22 24 26 27 4 5 18 19 20 21 25 12
## 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 4 4 5
## 13
## 5
summary(specnumber(mc_todas$Marzo)) # Resumen estadístico
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.000 2.000 2.000 2.333 3.000 5.000
#' ### Riqueza numérica de toda la "comunidad"
specnumber(colSums(mc_todas$Marzo))
## [1] 12
# TBI
mc_marzo_para_tbi <- cbind(
rbind(rownames(mc_pooled)[!rownames(mc_pooled) %in% rownames(mc_marzo)], mc_marzo) %>% sapply(., as.integer),
mc_mayo[,!colnames(mc_mayo) %in% colnames(mc_marzo)] %>% mutate_all(.funs = ~ifelse(.>0, 0, 0)))
colnames(mc_marzo)[!colnames(mc_marzo) %in% colnames(mc_mayo)] #Presentes en marzo, ausentes en mayo
## [1] "Pheidole moerens" "Tetramorium bicarinatum"
colnames(mc_mayo)[!colnames(mc_mayo) %in% colnames(mc_marzo)] #Presentes en mayo, ausentes en marzo
## [1] "Brachymyrmex heeri" "Pogonomyrmex schmitti" "Monomorium floricola"
## [4] "Acropyga dubitata"
mc_mayo_para_tbi <- cbind(
mc_mayo,
mc_marzo_para_tbi[,!colnames(mc_marzo_para_tbi) %in% colnames(mc_mayo)] %>% mutate_all(.funs = ~ifelse(.>0, 0, 0)))
mc_mayo_para_tbi <- mc_mayo_para_tbi[,colnames(mc_marzo_para_tbi)]
all.equal(colnames(mc_marzo_para_tbi), colnames(mc_mayo_para_tbi))
## [1] TRUE
rownames(mc_marzo_para_tbi) == rownames(mc_mayo_para_tbi)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
tbi_porc_diferencias <- TBI(mc_marzo_para_tbi, mc_mayo_para_tbi, method = '%difference', nperm = 999, test.t.perm = T)
tbi_porc_diferencias
## $TBI
## [1] 0.7142857 0.5000000 1.0000000 0.3333333 0.3333333 0.5555556 1.0000000
## [8] 0.5000000 0.6000000 0.3333333 0.5000000 0.4285714 0.3333333 1.0000000
## [15] 0.5000000 0.5000000 0.0000000 0.2000000 0.6000000 0.3333333 0.3333333
## [22] 0.2000000 1.0000000 0.6000000 0.4285714 0.0000000 0.3333333 0.6000000
##
## $p.TBI
## [1] 0.244 0.592 0.153 0.820 0.803 0.428 0.175 0.576 0.404 0.828 0.574 0.671
## [13] 0.820 0.171 0.568 0.538 0.997 0.912 0.412 0.828 0.823 0.931 0.173 0.382
## [25] 0.645 1.000 0.805 0.397
##
## $p.adj
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## $BCD.mat
## B/(2A+B+C) C/(2A+B+C) D=(B+C)/(2A+B+C) Change
## Site.1 0.7142857 0.0000000 0.7142857 -
## Site.2 0.2500000 0.2500000 0.5000000 0
## Site.3 0.5000000 0.5000000 1.0000000 0
## Site.4 0.1666667 0.1666667 0.3333333 0
## Site.5 0.1666667 0.1666667 0.3333333 0
## Site.6 0.0000000 0.5555556 0.5555556 +
## Site.7 0.2222222 0.7777778 1.0000000 +
## Site.8 0.2500000 0.2500000 0.5000000 0
## Site.9 0.2000000 0.4000000 0.6000000 +
## Site.10 0.0000000 0.3333333 0.3333333 +
## Site.11 0.0000000 0.5000000 0.5000000 +
## Site.12 0.4285714 0.0000000 0.4285714 -
## Site.13 0.2222222 0.1111111 0.3333333 -
## Site.14 0.5000000 0.5000000 1.0000000 0
## Site.15 0.0000000 0.5000000 0.5000000 +
## Site.16 0.2500000 0.2500000 0.5000000 0
## Site.17 0.0000000 0.0000000 0.0000000 0
## Site.18 0.2000000 0.0000000 0.2000000 -
## Site.19 0.4000000 0.2000000 0.6000000 -
## Site.20 0.1666667 0.1666667 0.3333333 0
## Site.21 0.1111111 0.2222222 0.3333333 +
## Site.22 0.0000000 0.2000000 0.2000000 +
## Site.23 0.5000000 0.5000000 1.0000000 0
## Site.24 0.2000000 0.4000000 0.6000000 +
## Site.25 0.2857143 0.1428571 0.4285714 -
## Site.26 0.0000000 0.0000000 0.0000000 0
## Site.27 0.0000000 0.3333333 0.3333333 +
## Site.28 0.0000000 0.6000000 0.6000000 +
##
## $BCD.summary
## mean(B/den) mean(C/den) mean(D) B/(B+C) C/(B+C) Change
## 0.2047902 0.2866497 0.4914399 0.4167147 0.5832853 +
##
## $t.test_B.C
## mean(C-B) Stat p.param p.perm p<=0.05
## Paired t.test 0.08185941 1.437708 0.1620069 0.174
##
## $BC
## [1] NA
##
## attr(,"class")
## [1] "TBI"
tbi_jaccard <- TBI(mc_marzo_para_tbi, mc_mayo_para_tbi, method = 'jaccard', nperm = 999, test.t.perm = T)
tbi_jaccard
## $TBI
## [1] 0.8333333 0.6666667 1.0000000 0.5000000 0.5000000 0.7142857 1.0000000
## [8] 0.6666667 0.7500000 0.5000000 0.6666667 0.6000000 0.5000000 1.0000000
## [15] 0.6666667 0.6666667 0.0000000 0.3333333 0.7500000 0.5000000 0.5000000
## [22] 0.3333333 1.0000000 0.7500000 0.6000000 0.0000000 0.5000000 0.7500000
##
## $p.TBI
## [1] 0.210 0.582 0.171 0.818 0.806 0.395 0.166 0.553 0.397 0.820 0.580 0.656
## [13] 0.795 0.170 0.556 0.554 0.999 0.919 0.435 0.817 0.829 0.923 0.163 0.398
## [25] 0.644 0.999 0.825 0.400
##
## $p.adj
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## $BCD.mat
## B/(A+B+C) C/(A+B+C) D=(B+C)/(A+B+C) Change
## Site.1 0.8333333 0.0000000 0.8333333 -
## Site.2 0.3333333 0.3333333 0.6666667 0
## Site.3 0.5000000 0.5000000 1.0000000 0
## Site.4 0.2500000 0.2500000 0.5000000 0
## Site.5 0.2500000 0.2500000 0.5000000 0
## Site.6 0.0000000 0.7142857 0.7142857 +
## Site.7 0.2222222 0.7777778 1.0000000 +
## Site.8 0.3333333 0.3333333 0.6666667 0
## Site.9 0.2500000 0.5000000 0.7500000 +
## Site.10 0.0000000 0.5000000 0.5000000 +
## Site.11 0.0000000 0.6666667 0.6666667 +
## Site.12 0.6000000 0.0000000 0.6000000 -
## Site.13 0.3333333 0.1666667 0.5000000 -
## Site.14 0.5000000 0.5000000 1.0000000 0
## Site.15 0.0000000 0.6666667 0.6666667 +
## Site.16 0.3333333 0.3333333 0.6666667 0
## Site.17 0.0000000 0.0000000 0.0000000 0
## Site.18 0.3333333 0.0000000 0.3333333 -
## Site.19 0.5000000 0.2500000 0.7500000 -
## Site.20 0.2500000 0.2500000 0.5000000 0
## Site.21 0.1666667 0.3333333 0.5000000 +
## Site.22 0.0000000 0.3333333 0.3333333 +
## Site.23 0.5000000 0.5000000 1.0000000 0
## Site.24 0.2500000 0.5000000 0.7500000 +
## Site.25 0.4000000 0.2000000 0.6000000 -
## Site.26 0.0000000 0.0000000 0.0000000 0
## Site.27 0.0000000 0.5000000 0.5000000 +
## Site.28 0.0000000 0.7500000 0.7500000 +
##
## $BCD.summary
## mean(B/den) mean(C/den) mean(D) B/(B+C) C/(B+C) Change
## 0.2549603 0.3610261 0.6159864 0.4139058 0.5860942 +
##
## $t.test_B.C
## mean(C-B) Stat p.param p.perm p<=0.05
## Paired t.test 0.1060658 1.455444 0.1570758 0.152
##
## $BC
## [1] NA
##
## attr(,"class")
## [1] "TBI"
tbi_sorensen <- TBI(mc_marzo_para_tbi, mc_mayo_para_tbi, method = 'sorensen', nperm = 999, test.t.perm = T)
tbi_sorensen
## $TBI
## [1] 0.7142857 0.5000000 1.0000000 0.3333333 0.3333333 0.5555556 1.0000000
## [8] 0.5000000 0.6000000 0.3333333 0.5000000 0.4285714 0.3333333 1.0000000
## [15] 0.5000000 0.5000000 0.0000000 0.2000000 0.6000000 0.3333333 0.3333333
## [22] 0.2000000 1.0000000 0.6000000 0.4285714 0.0000000 0.3333333 0.6000000
##
## $p.TBI
## [1] 0.219 0.580 0.170 0.805 0.814 0.420 0.157 0.564 0.415 0.819 0.554 0.643
## [13] 0.814 0.167 0.560 0.565 1.000 0.940 0.435 0.827 0.830 0.923 0.148 0.410
## [25] 0.622 1.000 0.789 0.400
##
## $p.adj
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## $BCD.mat
## B/(2A+B+C) C/(2A+B+C) D=(B+C)/(2A+B+C) Change
## Site.1 0.7142857 0.0000000 0.7142857 -
## Site.2 0.2500000 0.2500000 0.5000000 0
## Site.3 0.5000000 0.5000000 1.0000000 0
## Site.4 0.1666667 0.1666667 0.3333333 0
## Site.5 0.1666667 0.1666667 0.3333333 0
## Site.6 0.0000000 0.5555556 0.5555556 +
## Site.7 0.2222222 0.7777778 1.0000000 +
## Site.8 0.2500000 0.2500000 0.5000000 0
## Site.9 0.2000000 0.4000000 0.6000000 +
## Site.10 0.0000000 0.3333333 0.3333333 +
## Site.11 0.0000000 0.5000000 0.5000000 +
## Site.12 0.4285714 0.0000000 0.4285714 -
## Site.13 0.2222222 0.1111111 0.3333333 -
## Site.14 0.5000000 0.5000000 1.0000000 0
## Site.15 0.0000000 0.5000000 0.5000000 +
## Site.16 0.2500000 0.2500000 0.5000000 0
## Site.17 0.0000000 0.0000000 0.0000000 0
## Site.18 0.2000000 0.0000000 0.2000000 -
## Site.19 0.4000000 0.2000000 0.6000000 -
## Site.20 0.1666667 0.1666667 0.3333333 0
## Site.21 0.1111111 0.2222222 0.3333333 +
## Site.22 0.0000000 0.2000000 0.2000000 +
## Site.23 0.5000000 0.5000000 1.0000000 0
## Site.24 0.2000000 0.4000000 0.6000000 +
## Site.25 0.2857143 0.1428571 0.4285714 -
## Site.26 0.0000000 0.0000000 0.0000000 0
## Site.27 0.0000000 0.3333333 0.3333333 +
## Site.28 0.0000000 0.6000000 0.6000000 +
##
## $BCD.summary
## mean(B/den) mean(C/den) mean(D) B/(B+C) C/(B+C) Change
## 0.2047902 0.2866497 0.4914399 0.4167147 0.5832853 +
##
## $t.test_B.C
## mean(C-B) Stat p.param p.perm p<=0.05
## Paired t.test 0.08185941 1.437708 0.1620069 0.151
##
## $BC
## [1] NA
##
## attr(,"class")
## [1] "TBI"
# Calculando Jaccard de forma independiente
betadiver(rbind(colSums(mc_marzo_para_tbi), colSums(mc_mayo_para_tbi)), method = 'j')
## 1
## 2 0.625
# Para contrastar con tabla de valores significativos.
Resumen de resultado sobre TBI y diversidad Beta estacional
Calculando distintos indices de disimilaridad (%diferencias, Jaccard, Sorensen), y aplicando pruebas de permutacion a la diversidad Beta temporal (T1 y T2), se obtiene un resultado consistente sobre un cambio positivo en la diversidad de especies (ganancia de especies), pero este cambio no es significativo en terminos estadisticos.
Asimismo, tras calcular el indice de Jaccard mediante betadiver, contrastamos con la tabla de valores significativos de Raymundo Real, evidenciamos que existen diferencias significativas en la diversidad Beta comparando las matrices combinadas (pooled) de las estaciones para un nivel de significancia de 0.05, pero no para un nivel de significancia de 0.01. Esto sugiere que las diferencias en cuanto a la composicion de hormigas nidificantes entre ambas estaciones son realmente bajas.
precip_orig <- read_excel('JARDIN BOTANICO (JOSE RAMON).xls', range = "A9:I40")
precip <- precip_orig %>%
mutate_all(~ gsub('INAP', '0', .x)) %>%
mutate_all(~ gsub('-', NA, .x)) %>%
mutate_all(~ as.numeric(.x))
sapply(precip, summary)
## $DIA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.0 8.5 16.0 16.0 23.5 31.0
##
## $ENE
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2484 0.0000 3.3000
##
## $FEB
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 0.000 0.000 1.518 1.925 11.400 3
##
## $MAR
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 1.000 7.116 13.200 34.800
##
## $ABR
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 0.000 0.000 7.323 6.025 63.200 1
##
## $MAY
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 3.035 2.250 29.500
##
## $JUN
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 0.000 0.000 3.843 1.300 68.800 1
##
## $JUL
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 1.000 9.097 13.950 47.800
##
## $AGO
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 0.000 0.000 7.353 17.325 34.500 1
sapply(precip, sum, na.rm=T)
## DIA ENE FEB MAR ABR MAY JUN JUL AGO
## 496.0 7.7 42.5 220.6 219.7 94.1 115.3 282.0 220.6
precip_mar_colec <- precip %>%
select(DIA, MAR) %>%
filter(DIA %in% c(21:25, 30, 31)) %>%
pull(MAR)
sum(precip_mar_colec)
## [1] 27.1
precip_may_colec <- precip %>%
select(DIA, MAY) %>%
filter(DIA %in% c(10:13, 16:19, 31)) %>%
pull(MAY)
sum(precip_may_colec)
## [1] 46.7
library(sp)
# Matriz transformada Hellinger
mc_nidos_hell <- decostand(mc_nidos, 'hellinger')
# Matriz ambiental con centroides
st_geometry(hex_inters) <- 'geom'
ma_sf <- hex_inters %>%
inner_join(ma %>%
rownames_to_column('ENLACE') %>%
mutate(ENLACE = as.integer(ENLACE)))
ma_sp <- ma_sf %>% as_Spatial()
centroides <- ma_sf %>% st_centroid
ma_xy <- centroides %>% st_coordinates %>% as.data.frame
(vecindad <- centroides %>% dnearneigh(d1 = 0, d2 = 75))
## Neighbour list object:
## Number of regions: 28
## Number of nonzero links: 128
## Percentage nonzero weights: 16.32653
## Average number of links: 4.571429
plot(ma_sp)
plot(vecindad, coords = ma_xy, add=T, col = 'red')
(pesos_b <- nb2listw(vecindad, style = 'B'))
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 28
## Number of nonzero links: 128
## Percentage nonzero weights: 16.32653
## Average number of links: 4.571429
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 28 784 128 256 2544
# Matriz sin tendencia
mc_nidos_sin_tendencia <- resid(
lm(as.matrix(mc_nidos_hell) ~ .,
data = ma_xy))
(autocor_global_residuos_mc <- sapply(
dimnames(mc_nidos_sin_tendencia)[[2]],
function(x)
moran.mc(
x = mc_nidos_sin_tendencia[,x],
listw = pesos_b,
zero.policy = T,
nsim = 9999),
simplify = F))
## $`Pheidole subarmata`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = -0.12312, observed rank = 2306, p-value = 0.7694
## alternative hypothesis: greater
##
##
## $`Dorymyrmex antillanus`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = -0.084899, observed rank = 3532, p-value = 0.6468
## alternative hypothesis: greater
##
##
## $`Paratrechina longicornis`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = -0.13807, observed rank = 1794, p-value = 0.8206
## alternative hypothesis: greater
##
##
## $`Solenopsis geminata`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = 0.14283, observed rank = 9320, p-value = 0.068
## alternative hypothesis: greater
##
##
## $`Pheidole jelskii`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = -0.16881, observed rank = 1021, p-value = 0.8979
## alternative hypothesis: greater
##
##
## $`Pheidole moerens`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = -0.093401, observed rank = 1420, p-value = 0.858
## alternative hypothesis: greater
##
##
## $`Mycetomoellerius jamaicensis`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = -0.082011, observed rank = 3406, p-value = 0.6594
## alternative hypothesis: greater
##
##
## $`Tetramorium lanuginosum`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = -0.093645, observed rank = 2933, p-value = 0.7067
## alternative hypothesis: greater
##
##
## $`Odontomachus ruginodis`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = -0.29487, observed rank = 30, p-value = 0.997
## alternative hypothesis: greater
##
##
## $`Wasmannia auropunctata`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = -0.093884, observed rank = 3161, p-value = 0.6839
## alternative hypothesis: greater
##
##
## $`Cyphomyrmex rimosus`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = 0.024807, observed rank = 7402, p-value = 0.2598
## alternative hypothesis: greater
##
##
## $`Tetramorium bicarinatum`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = -0.097105, observed rank = 1337, p-value = 0.8663
## alternative hypothesis: greater
##
##
## $`Brachymyrmex heeri`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = -0.096045, observed rank = 287, p-value = 0.9713
## alternative hypothesis: greater
##
##
## $`Pogonomyrmex schmitti`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = -0.19293, observed rank = 531, p-value = 0.9469
## alternative hypothesis: greater
##
##
## $`Monomorium floricola`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = 0.018724, observed rank = 8068, p-value = 0.1932
## alternative hypothesis: greater
##
##
## $`Acropyga dubitata`
##
## Monte-Carlo simulation of Moran I
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
## number of simulations + 1: 10000
##
## statistic = 0.018724, observed rank = 8091, p-value = 0.1909
## alternative hypothesis: greater
(autocor_global_residuos <- sapply(
dimnames(mc_nidos_sin_tendencia)[[2]],
function(x)
moran.test(
x = mc_nidos_sin_tendencia[,x],
listw = pesos_b,
zero.policy = T),
simplify = F))
## $`Pheidole subarmata`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = -0.77838, p-value = 0.7818
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.12312262 -0.03703704 0.01223152
##
##
## $`Dorymyrmex antillanus`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = -0.41818, p-value = 0.6621
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.08489939 -0.03703704 0.01309983
##
##
## $`Paratrechina longicornis`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = -0.92421, p-value = 0.8223
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.13807459 -0.03703704 0.01195157
##
##
## $`Solenopsis geminata`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = 1.576, p-value = 0.05752
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.14282977 -0.03703704 0.01302614
##
##
## $`Pheidole jelskii`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = -1.1842, p-value = 0.8818
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.16880693 -0.03703704 0.01238092
##
##
## $`Pheidole moerens`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = -1.0731, p-value = 0.8584
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.093401039 -0.037037037 0.002759073
##
##
## $`Mycetomoellerius jamaicensis`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = -0.50945, p-value = 0.6948
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.082011111 -0.037037037 0.007793372
##
##
## $`Tetramorium lanuginosum`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = -0.62151, p-value = 0.7329
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.093645079 -0.037037037 0.008295785
##
##
## $`Odontomachus ruginodis`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = -2.3369, p-value = 0.9903
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.29487015 -0.03703704 0.01217280
##
##
## $`Wasmannia auropunctata`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = -0.50626, p-value = 0.6937
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.09388384 -0.03703704 0.01260874
##
##
## $`Cyphomyrmex rimosus`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = 0.56388, p-value = 0.2864
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.02480714 -0.03703704 0.01202887
##
##
## $`Tetramorium bicarinatum`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = -1.1375, p-value = 0.8723
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.097105008 -0.037037037 0.002788338
##
##
## $`Brachymyrmex heeri`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = -1.904, p-value = 0.9715
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.0960451331 -0.0370370370 0.0009604943
##
##
## $`Pogonomyrmex schmitti`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = -1.4674, p-value = 0.9289
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.19292906 -0.03703704 0.01128672
##
##
## $`Monomorium floricola`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = 0.88658, p-value = 0.1877
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.018724330 -0.037037037 0.003955756
##
##
## $`Acropyga dubitata`
##
## Moran I test under randomisation
##
## data: mc_nidos_sin_tendencia[, x]
## weights: pesos_b
##
## Moran I statistic standard deviate = 0.88658, p-value = 0.1877
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.018724330 -0.037037037 0.003955756
pvalor_moran <- sapply(names(autocor_global_residuos),
function(x) autocor_global_residuos[[x]]$p.value)
tabla_patrones_espaciales <- data.frame(
Especie = names(pvalor_moran),
ValorP = pvalor_moran,
Patron = ifelse(pvalor_moran<0.1, 'Agregado', 'Aleatorio'))
tabla_patrones_espaciales %>% kable
| Especie | ValorP | Patron | |
|---|---|---|---|
| Pheidole subarmata | Pheidole subarmata | 0.7818267 | Aleatorio |
| Dorymyrmex antillanus | Dorymyrmex antillanus | 0.6620916 | Aleatorio |
| Paratrechina longicornis | Paratrechina longicornis | 0.8223113 | Aleatorio |
| Solenopsis geminata | Solenopsis geminata | 0.0575183 | Agregado |
| Pheidole jelskii | Pheidole jelskii | 0.8818410 | Aleatorio |
| Pheidole moerens | Pheidole moerens | 0.8583758 | Aleatorio |
| Mycetomoellerius jamaicensis | Mycetomoellerius jamaicensis | 0.6947807 | Aleatorio |
| Tetramorium lanuginosum | Tetramorium lanuginosum | 0.7328685 | Aleatorio |
| Odontomachus ruginodis | Odontomachus ruginodis | 0.9902783 | Aleatorio |
| Wasmannia auropunctata | Wasmannia auropunctata | 0.6936615 | Aleatorio |
| Cyphomyrmex rimosus | Cyphomyrmex rimosus | 0.2864181 | Aleatorio |
| Tetramorium bicarinatum | Tetramorium bicarinatum | 0.8723456 | Aleatorio |
| Brachymyrmex heeri | Brachymyrmex heeri | 0.9715441 | Aleatorio |
| Pogonomyrmex schmitti | Pogonomyrmex schmitti | 0.9288624 | Aleatorio |
| Monomorium floricola | Monomorium floricola | 0.1876521 | Aleatorio |
| Acropyga dubitata | Acropyga dubitata | 0.1876521 | Aleatorio |
# mc_nidos_sin_tendencia_d <- dist(mc_nidos_sin_tendencia)
# (mc_nidos_correlograma <- mantel.correlog(
# mc_nidos_sin_tendencia_d,
# XY = ma_xy,
# nperm = 999))
# plot(mc_nidos_correlograma)